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Reynolds-averaged and hybrid Reynolds-averaged/large-eddy simulations have been ap- 
plied to a supersonic coaxial jet flow experiment. The experiment utilized either helium 
or argon as the inner jet nozzle fluid, and the outer jet nozzle fluid consisted of labo- 
ratory air. The inner and outer nozzles were designed and operated to produce nearly 
pressure- matched Mach 1.8 flow conditions at the jet exit. The purpose of the computa- 
tional effort was to assess the state-of-the-art for each modeling approach, and to use the 
hybrid Reynolds-averaged/large-eddy simulations to gather insight into the deficiencies 
of the Reynolds-averaged closure models. The Reynolds- averaged simulations displayed 
a strong sensitivity to choice of turbulent Schmidt number. The baseline value chosen 
for this parameter resulted in an over-prediction of the mixing layer spreading rate for 
the helium case, but the opposite trend was noted when argon was used as the injectant. 
A larger turbulent Schmidt number greatly improved the comparison of the results with 
measurements for the helium simulations, but variations in the Schmidt number did not 
improve the argon comparisons. The hybrid simulation results showed the same trends 
as the baseline Reynolds-averaged predictions. The primary reason conjectured for the 
discrepancy between the hybrid simulation results and the measurements centered around 
issues related to the transition from a Reynolds-averaged state to one with resolved tur- 
bulent content. Improvements to the inflow conditions are suggested as a remedy to this 
dilemma. Comparisons between resolved second-order turbulence statistics and their mod- 
eled Reynolds-averaged counterparts were also performed. 
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Cartesian velocity components 
Cartesian coordinates 
cylindrical coordinates 

non-dimensional “law of the wall” coordinate 

mass fraction of species m 

blending function constant 

Kronecker delta 

sub grid scale filter width 

blending function length scale ratio 

von Karman constant or MUSCL parameter 

molecular viscosity 

turbulent viscosity 

kinematic viscosity 

density 

specific turbulence dissipation rate 


Subscripts 

i,j,k computational coordinate indices 

inj center jet inject ant index 

air coflow air index 


Introduction 

Reynolds- averaged Computational Fluid Dynamics (CFD) models have become an integral part of the 
design and analysis of high-speed air-breathing engines. The maturation of multi-purpose CFD codes coupled 
with advancements in computer architectures have substantially reduced the turn-around time required to 
perform steady-state Reynolds- Averaged Simulations (RAS). Unfortunately, the turbulence models required 
to close the Reynolds-averaged equation set have not kept pace with these advancements. As a result, RAS 
approaches often require calibrations for key model parameters. One example is given in Ref. 1, where 
variations of the turbulent Prandtl and Schmidt numbers for a scramjet combustor simulation were shown 
to produce outcomes that ranged from engine unstart to complete flame blow-out. Similar examples that 
illustrate solution sensitivities to unknown RAS modeling parameters can be found in Refs. 2 and 3. 

Large Eddy Simulation (LES) methods have the potential to reduce the modeling sensitivity inherent to 
RAS approaches, since the intent of LES is to resolve the large scale turbulent structures while modeling 
only the small scales. Regrettably, the computational expense of this approach (particularly when applied 
to configurations of interest to the high-speed propulsion community) is well beyond what can be deemed 
as practical by today’s standards. Hybrid RAS/LES approaches offer some relief to the computational costs 
associated with LES. These methodologies allow LES content to be resolved in areas that require a more 
rigorous modeling approach, while maintaining a more cost effective RAS approach for benign regions of the 
flow (e.g. attached boundary layers). Hybrid approaches first appeared in the literature a decade ago, 4 ’ 5 
and while many advancements have been made since then, the modeling is far from mature as evident by 
the vast array of hybrid approaches that have appeared in recent years. 6 ’ 7 ’ 8? 9 The computational expense 
required for a hybrid simulation, while less than that of a full LES, is still formidable when compared with 
steady- state RAS. 

The present effort utilizes a hybrid RAS/LES approach to model a series of supersonic coaxial jet ex- 
periments 10, 11 that have been studied at the NASA Langley Research Center. These experiments involve 
the coaxial injection of either helium or argon into air, and CFD data (based on RAS approaches) were also 
obtained and compared to measurements in the works referenced above. The goals of this computational 
effort are to: 

• Assess the capabilities of RAS and hybrid RAS/LES for predicting high speed mixing layers 

• Determine the model sensitivities for both RAS and hybrid RAS/LES approaches 

• Attempt to use the hybrid RAS/LES results to assess the appropriateness of models used by RAS 

Measured values of composition, Pitot pressure, velocity, and rms of the velocity fluctuation are available 
for comparison with the simulations. Additional higher-order correlations (that are difficult to obtain exper- 
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imentally) are extracted from the hybrid RAS/LES results and compared with the RAS predictions. The 
correlations considered include the Reynolds stress tensor and the Reynolds mass flux vector. 

Geometry Description and Flow Conditions 

A schematic of the coaxial nozzle assembly is shown in Fig. 1. The injectant supplied to the center jet 
nozzle was either a mixture of 95% helium and 5% oxygen (by volume), or pure argon. In the former case, 
a trace amount of oxygen was added to allow for the measurement of the streamwise component of velocity 
using the RELIEF 12 oxygen flow-tagging technique. The internal diameter of the center jet nozzle is 10 mm 
at the nozzle exit. The center-body that forms the internal nozzle is 0.25 mm thick at the nozzle exit, 
providing a small blunt base to anchor the shear layer formed between the two nozzle streams. The coflow 
air nozzle has an internal diameter of 60.47 mm, and the outer surface of this nozzle extends 12.66 mm 
farther downstream than the center-body. The 38.6° juncture between the internal surface of the coflow 
nozzle and the conical exterior surface is sharp. Hence, there is no appreciable base region to segregate the 
outer jet flow from the surrounding ambient air. Further details concerning the geometry of the rig, and the 
methodology used for its design, can be found in Ref. 10. 



AH dimensions in mm 


Figure 1: Schematic of the coaxial nozzle assembly 

The nominal flow conditions for each experiment are given in Tables 1 and 2 along with the cited 
uncertainties. 10 ’ 11 The thermocouples used to measure the jet temperatures were located in the gas supply 
lines, and the pressure taps used to measure the operating pressures were positioned as shown in Fig. 1. 
Although each jet has a design Mach number of 1.8, the flow velocities are markedly different due to the 
difference in the molecular weights for each center jet injectant. The velocity of the center jet is more than 
twice that of the coflow jet for the helium case (Case 1), but it is 16 % lower than that of the coflow jet for 
the argon case (Case 2). An estimation of the convective Mach number: 

'Uair 'U j inj 

C^air T ttinj 

yields a value of 0.7 for Case 1 and 0.16 for Case 2, respectively. Hence, one expects compressibility effects 
to be evident in Case 1, while Case 2 is expected to behave more like an incompressible shear layer. 
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Table 1: Case 1 - Helium- Air Test Conditions 


Nominal Conditions 

Center Jet 

Coflow Jet 

Ambient 

Mach Number 

1.8 a 

1.8 a 

0.025 6 

Total Temperature [K] 

305.0 (±15) 

300.0 (±6) 

294.6 (±6) 

Total Pressure [kPa] 

614.93 (±6) 

579.80 (±4) 

101.325 (±1) 


a nozzle design Mach number, b value assumed for the entrained ambient flow 


Table 2: Case 2 - Argon-Air Test Conditions 


Nominal Conditions 

Center Jet 

Coflow Jet 

Ambient 

Mach Number 

1.8 a 

1.8 a 

0.025 b 

Total Temperature [K] 

297.9 (±3.5) 

294.3 (±3.5) 

294.6 (±3.5) 

Total Pressure [kPa] 

615.86 (±5.5) 

580.68 (±4.4) 

101.325 (±0.6) 


a nozzle design Mach number, b value assumed for the entrained ambient flow 


Computational Methodology 

All computational results were obtained using the VULCAN (Viscous Upwind aLgorithm for Complex 
flow ANalysis) software package. The code solves the unsteady, conservation equations appropriate for 
calorically or thermally perfect gases with a cell-centered finite volume scheme. Efficient utilization of 
parallel architectures is realized through calls to MPI (Message Passing Interface) routines using an SPMD 
(Single Program Multiple Data) paradigm; a natural choice for multi-block flow solvers. Arbitrary block- 
to-block (non-aligned) connectivity is also available, allowing the flexibility to add or remove grid points 
at zonal interfaces. A variety of upwind formulations are available for evaluating the inviscid fluxes, while 
central differences are used for the viscous fluxes. Diagonalized Approximate Factorization (DAF) or planar 
relaxation (ILU) are the primary options available for steady- state simulations, while an explicit multi-stage 
Runge-Kutta or an implicit dual time-stepping strategy (utilizing DAF or ILU) are the options available for 
unsteady applications. A variety of one-equation and two-equation turbulence models exist for Reynolds- 
averaged simulations. LES and hybrid RAS/LES models are also available for turbulent flows that require 
a more rigorous modeling effort. Further details describing the code are found elsewhere. 13, 14 

Physical and Numerical Model Description 

Both RAS and hybrid RAS/LES computational approaches were considered in this effort. The RAS 
results utilized the Low-Diffusion Flux Split Scheme (LDFSS) of Edwards 15 to evaluate the inviscid fluxes. 
The Monotone Upstream-centered Scheme for Conservation Laws (MUSCL) extrapolation parameter (k) 
was chosen as 1/3 to minimize spatial truncation error, and the van Leer flux limiter 16 was employed 
to enforce the Total Variation Diminishing (TVD) property. All of the steady- state RAS solutions were 
advanced in time with the diagonalized approximate factorization scheme. The turbulence model chosen 
was the Wilcox (1998) k-cu model, 17 and the wall function procedure of Wilcox 18 was used to relax the 
grid spacing requirements near solid surfaces. A dilatation-dissipation modification 17 to the turbulence 
kinetic energy equation was also enabled to reduce the shear layer spreading rate that has been observed 
at high compressible Mach number conditions. The baseline values for the turbulent Prandtl ( Pr t ) and 
Schmidt (Sc t ) numbers, which control the turbulent transport of energy and mass, were chosen as 0.9 and 
0.5, respectively. The turbulent Schmidt number was one of the properties that was varied in parametric 
studies to evaluate the sensitivity of the RAS solutions to the value assumed for this coefficient. 

The hybrid RAS/LES results required a somewhat different algorithmic approach to encourage the de- 
velopment of resolved turbulent content. In particular, the standard Total Variation Diminishing (TVD) 
limiters that are typically employed for RAS are overly dissipative when used for LES. This class of lim- 
iter tends to locally reduce the accuracy of the inviscid flux scheme to first-order whenever extrema of any 
magnitude are encountered. Essentially Non-Oscillatory (ENO) limiters offer a simple means of maintaining 
second-order accuracy near local extrema. This is accomplished by utilizing data at additional grid nodes to 
alter the stencil used to construct the inviscid flux terms. The ENO limiter used in this effort is the SONIC- 
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A limiter of Suresh and Huynh, 19 which can be considered as a second-order extension of the van Leer 
TVD limiter. All of the time-accurate hybrid RAS/LES solutions were advanced in time using a dual time- 
stepping approach that combined the DAF method for integration in pseudo-time, with a 3-point backwards 
finite difference approximation for integration in real-time. The values selected for the physical time-step 
and sub-iteration CFL constraint were 0.05 {is and 50.0, respectively. The time-step was chosen based on a 
cell residence time constraint to ensure that turbulent structures would not traverse more than one grid cell 
per iteration. The sub-iteration process was considered converged when the residual error dropped at least 
2.5 orders of magnitude. This level of convergence typically required 4-7 sub-iterations for each physical 
time- step. 

The hybrid RAS/LES methodology used in this effort builds on the previous work discussed in Ref. 6. 
This framework was designed to enforce a RAS behavior near solid surfaces, and to switch to an LES 
behavior in the outer portion of the boundary layer and free shear regions. Hence, this formulation can be 
thought of as a wall-modeled LES approach, where RAS is used as the near-wall model. The basic idea is 
to blend any trusted RAS eddy viscosity with a desired LES Sub-Grid Scale (SGS) viscosity, along with any 
transport equations that involve a common RAS and SGS property. In this effort, the Wilcox (1998) k-u 
RAS model 17 was blended with the one-equation model of Yoshizawa. 20 The Yoshizawa model involves an 
evolution equation for the SGS turbulence kinetic energy, hence the blended expressions that are appropriate 
for this model combination are: 

Hybrid RAS/SGS viscosity = (F) [RAS viscosity] + (1 — F) [SGS viscosity] 

Hybrid RAS/SGS k-equation = (F) [RAS k-equation] + (1 — F) [SGS k-equation] (2) 

where F is a blending function that varies between 0 and 1. Note that the transport equation for the RAS 
specific dissipation rate (uj) does not have an SGS counterpart. Hence, the blending is not applied to this 
equation, and all of the terms in this equation that involve the eddy viscosity are evaluated based on RAS 
relationships. 

The motivation behind the development of this hybrid RAS/LES framework is two- fold. First, the 
blending of two independent RAS and LES closure models offers the flexibility of having an optimized set of 
closure equations for both RAS and LES modes. The second (and more critical) driving factor was the desire 
to alleviate the difficulties associated with the design of grid topologies that are appropriate for purely grid- 
dependent blending paradigms such as those used for Detached Eddy Simulation (DES). 4 ’ 21 The movement 
away from simple grid-dependent blending strategies has gained momentum in recent years. 7 ’ 8 In fact, even 
the developers of DES are now promoting an “improved” version of their scheme termed Delayed Detached 
Eddy Simulation (DDES) 8 that involves flow-dependent functions to switch between RAS and LES regimes. 

The blending function (F) used in this effort is based on the ratio of the wall distance d to a modeled 
form of the Taylor microscale (£): 



where k is the von Karman constant (0.41), C M is 0.09, a is a user-defined model constant, and is set to 
tan -1 (0.98) to force the balancing position of F (i.e. the position where k rj 2 = y / C^) to 0.99. The value 
chosen for a provides control over the y + position where the average LES to RAS transition point (defined 
as F = 0.99) occurs. If resolved LES content is desired for an attached boundary layer, then this constant 
should be set such that the transition point occurs in the region where the boundary layer wake law starts 
to deviate from the log law. If the transition point is enforced at a lower y + value (e.g. well within the log 
law region), a dual log layer can appear. 22, 23 Conversely, if the transition point is enforced at a y + value 
that extends well into the defect layer, then the level of resolved turbulence will be limited. Details on a 
procedure to analytically determine the value for a that corresponds to a target y + value is described in 
Ref. 23. 
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Grid Details 

A two-dimensional (axisymmetric) grid was generated for the steady-state Reynolds-averaged solutions 
(see Fig. 2). This grid consisted of just under 250,000 cells divided across five structured grid zones. The 
computational domain extended 150 center jet diameters (D) downstream of the jet exit, and between 40 
and 70 jet diameters in the radial direction. The internal flowfield within the concentric nozzles was also 
included as part of the computational domain to allow for the development of boundary layers. The last 
experimental data stations were located 26.101 (Case 1) and 45.276 (Case 2) jet diameters downstream of 
the center jet exit plane. Hence, most of the grid was clustered between the jet exit and an x/D of 50. The 
portion of the computational domain downstream of x/D = 50 was added to provide plenty of space for 
the jet flow to reach a subsonic state, which ensured that the specified pressure outflow boundary condition 
remained well-posed. The far-held boundary was also placed many tens of jet diameters away from the 
domain of interest to minimize any chance of data corruption that might occur via wave reflections off of 
the characteristic inflow condition that was applied along this boundary. The nozzle stagnation conditions 
given in Tables 1 and 2 were applied at each nozzle inflow plane. The operating total temperature of each 
nozzle was nominally the same as the room temperature, so all solid surfaces were assumed to be adiabatic. 
The grid was clustered to all solid surfaces at a level appropriate for the use of wall functions (t/ + < 36). 
Finally, the entire grid system was broken up into a total of 179 grid blocks which yielded excellent load 
balance statistics for up to 64 processors. 

An azimuthal slice of the three-dimensional grid generated for the hybrid RAS/LES cases is shown 
in Fig. 3. The boundary conditions and the extent of the computational domain was identical to that 
used for the axisymmetric RAS cases. The gridding strategy, however, was altered to reflect the change 
in computational algorithm. LES requires the use of grid cells that are roughly isotropic to allow for the 
development and sustainment of resolved turbulent structures. Hence, the streamwise grid spacing had to 
be reduced from that utilized for the pure RAS. This requirement, coupled with the fact that the full 3-D 
flowfield must be solved, forced a concession to be made to keep the computational costs within reason. The 
compromise that was accepted was to provide a level of grid resolution capable of supporting LES only up 
to an x/D of 25 (highlighted in red). This domain captures all but the last experimental station for Case 1, 
and 13 of the 16 stations for Case 2. Non-aligned zonal interfaces (patches) were inserted at x/D stations 
of 25, 50, and 100 to systematically remove grid nodes in the radial direction. A patch interface was also 
inserted in the ambient air region to coarsen the streamwise spacing in the far-held where axial refinement 
is not required. Each patch interface is denoted by a solid black line in Fig. 3. In addition to substantially 
reducing the computational costs, the coarsening process performed at large x/D tends to dissipate the large 
scale vortical structures that propagate out of the LES domain of interest. This reduces the likelihood of 
an ill-posed boundary condition appearing at the outflow plane. An H-0 mesh topology was chosen for the 
cross-how planes to avoid the numerical difficulties associated with collapsed faces that would result from a 
classic polar topology. This feature of the grid is illustrated in the three-dimensional view of the LES domain 
shown in Fig. 4. Finally, the entire grid system was broken up into a total of 1669 grid blocks which resulted 
in excellent load balance statistics for up to 360 processors. 

Reynolds- Averaged Results 

Contours of helium mass fraction and Mach number extracted from the baseline RAS (, Sc t = 0.5) are 
shown in Fig. 5. The helium mass fractions have been normalized by the the center jet value (95% by volume 
corresponds to 70.39% by mass): 

Y H e^Y He / 0.7039 (4) 

for all of the results that follow to simplify the analysis and comparisons with the argon cases. The helium 
contours show that the potential core of the center jet persists for approximately 12 jet diameters, and the 
center-line helium mass fraction in the plume that forms downstream of this station decays to a value of 
0.115 at an x/D = 50. The Mach contours show that both nozzle streams are nearly pressure matched. 
Although not evident at this scale, a pair of counter-rotating separation bubbles are present behind the small 
blunt base of the center-body. The weak expansion fan that is seen on either side of the base region forms as 
the jets negotiate around this separated flow zone. A weak recompression shock then develops as the flow is 
forced to turn back on itself at the close-off point of the recirculation bubble. The conical shock that forms 
in the center jet flowfield steepens as it approaches the axis of symmetry, resulting in a Mach disk at the 
axis. As will be shown later, the total pressure loss that occurs across this shock wave is also evident in the 
measured Pitot pressure surveys. 
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Figure 4: Isometric visualization of the LES-resolved portion of the hybrid RAS/LES grid 
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Figure 5: Contours of normalized helium mass fraction and Mach number for the baseline RAS 

Measured helium mass fraction profiles are compared with the RAS results in Fig. 6. The measurements 
of injectant composition were made using a gas sampling probe as detailed in Ref. 10. In general, the 
baseline simulation over-predicted the rate of mixing between the coaxial streams. Hence, an additional 
simulation was performed with the Sc t value increased by a factor of 2 (to 0.25). The third result shown in 
this figure is a repeat of the baseline condition, but simulated on a coarse grid (factor of 2 coarsening in each 
coordinate direction). Comparisons are made at 8 axial stations downstream of the center jet exit plane. 
The results extracted from the baseline (Sc t = 0.5) simulation over-predicted the measured growth rate of 
the mixing layer for all but the first experimental station. It should be noted that the RAS model used 
in this effort included a dilatation-dissipation modification 17 to reduce the shear layer spreading rate for 
compressible mixing layers. This leaves the turbulent Schmidt number (which controls the rate of turbulent 
mass diffusion) as the remaining parameter that can be adjusted to reduce the mixing rate. The results 
obtained when doubling the Sc t value (which reduces the modeled turbulent diffusion term by a factor of 
2) agree well with the measurements. In fact, the results obtained with Sc t = 1.0 matched nearly every 
measured result to within the precision of the measurement device of ±1.0 — 1.5%. The effect of coarsening 
the grid resolution had a much smaller impact on the predicted results, indicating that the level of grid 
convergence is well within the uncertainty bounds of the physical models employed. 

Measured Pitot pressure profiles 10 are compared with the RAS results in Fig. 7. The large deficit in 
Pitot pressure near r/D = 0.5 is a result of viscous losses from the upstream boundary layers. The small 
deficit seen near the axis of symmetry at the first four axial stations is a result of the Mach disk formed by 
the coalescence of the conical shock front that was discussed previously (see Fig. 5). The measured value 
of this deficit is slightly larger than the predicted value. The Pitot pressure deficit is also under-predicted 
for the coflow shock structure that appears near r/D = 1.1 at the first axial station. In general, the best 
agreement with measurements is given by the Sc t = 1.0 data. The level of agreement seen between this 
simulation and the measured values is within the measurement uncertainty (based solely on the transducer 
error) of ±0.5% at nearly every measured point for x/D < 15. Downstream of this station, the Sc t = 1.0 
Pitot data under-predicts the measured values near the axis, and over-predicts them near the lower edge of 
the shear layer present between the coflow jet and the ambient surroundings. Once again, the differences 
noted between the fine and coarse grid solutions is much smaller than that associated with a modification 
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of the turbulent Schmidt number. 

Contours of argon mass fraction and Mach number extracted from the baseline RAS (Sc t = 0.5) are 
shown in Fig. 8. The argon results show a much longer potential core than that for the helium case. In 
fact, the potential core in the simulations extends beyond the x/D = 50 station shown in Fig. 8. The slower 
mixing rate between the two jet streams is a result of the reduced shear associated with Case 2. The velocity 
of the argon stream is only 16% lower than the coflow stream. This is contrasted with the helium case, where 
the velocity of the helium stream was more than twice that of the coflow stream. The Mach contours clearly 
show that the nozzle streams are nearly pressure matched, and the near-field flow structures are practically 
identical to those discussed previously for Case 1. The far- field shock/expansion pattern is quite different, 
however, due to the reduced spreading rate of the mixing layer. 

Measured argon mass fraction profiles 11 are compared with the RAS results in Fig. 9. The baseline argon 
simulation under-predicted the rate of mixing between the coaxial streams. This is in contrast to the helium 
case where the baseline modeling parameters over-predicted the growth rate of the mixing layer. Hence, 
an additional simulation was performed with the turbulent Schmidt number reduced by a factor of 2. A 
third simulation was performed with the baseline modeling parameters using the coarse grid to evaluate the 
grid dependence. Comparisons are made at the same 8 axial stations used to evaluate the helium injection 
data. The results obtained when reducing the Sc t value had the desired effect of reducing the length of the 
potential core of the argon jet. However, this adjustment also resulted in a peculiar inflection point near 
the outer edge of the mixing layer that is not evident in the measurements. The reduced Schmidt number 
results also over-predicted the spreading rate of the mixing layer for x/D < 18. Overall, the best result was 
obtained with the baseline modeling parameters, which predicted the correct spreading rate of the shear layer 
for x/D < 18. The convective Mach number for the mixing layer was estimated to be 0.16, so an additional 
simulation (not shown) was performed with the dilatation-dissipation modification disabled to ensure that 
the model was not impacting the predictions. As expected, the compressibility correction had no impact 
on this simulation. One interesting feature noticed in the experimental data is an increased deviation from 
symmetry that begins to appear at the x/D = 12 station. This happens to be the station where simulations 
begin to deviate from the measurements. The experimental probes were traversed over the entire coaxial 
flowfield (i.e. data was gathered for both positive and negative y - values relative to the jet axis). Since the 
data should be axisymmetric in the mean, the measured data was plotted against the absolute value of y in 
this work, which effectively results in two sets of “radial” data at each experimental station. The deviation 
from symmetry is much larger for Case 2 than for Case 1, which might have implied a slight misalignment of 
the jet assembly during the argon tests. 11 As was seen for the helium injection conditions, the predictions 
offered by both the fine and coarse meshes were quite similar. 

Measured Pitot pressure profiles 11 for Case 2 are compared with the RAS results in Fig. 10. The near- 
field Pitot profile features are identical to those seen for Case 1, as implied by the Mach contours that were 
discussed earlier. In general, the level of agreement seen between the simulations and the measured values 
is within the measurement uncertainty for x/D < 15. Downstream of this station, the computed Pitot data 
shows sharper features, which is a direct result of the under-prediction of the mixing processes at this point 
in the flowfield. It is interesting to note that the over-prediction of the Pitot pressure near the lower edge 
of the shear layer present between the coflow jet and the ambient surroundings is larger that what was seen 
for Case 1. The shear between the coflow jet and the ambient environment is practically identical to that 
for the helium condition. Hence, the under-prediction of the argon/air mixing layer has evidently modified 
the shock system in the outer jet enough to alter the behavior of the coflow /ambient air shear layer. 


Hybrid RAS/LES Results 


The coaxial flowfields examined in this effort have three distinct flow regions: center jet, coflow jet, and 
the ambient surroundings. The core velocity in the coflow jet was 480 m/s, and the center jet values were 
1100 m/s (helium) and 400 m/s (argon). The ambient air region was almost stagnant. This presented a 
problem in that the hybrid RAS/LES approach requires a time-accurate integration with a time-step small 
enough to temporally resolve turbulent structures at scales representative of a computational cell width. A 
time-step appropriate for resolving a turbulent structure in the jet flow is given by the following estimate: 


At 


Ax 

max ( Uinj 5 ^air) 


(5) 
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Figure 6: Comparison of normalized helium mass fraction (RAS predictions) with measured values 
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Figure 7: Comparison of Case 1 Pitot pressure (RAS predictions) with measured values 
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Figure 8: Contours of argon mass fraction and Mach number for the baseline RAS 


Based on the disparate characteristic velocities between the jet flows and the ambient air, this time-step 
would require hundreds of iterations to transport a particle across one cell length in the ambient air region. 
Fortunately, the shear layer of interest is the mixing layer between the supersonic jets rather than the shear 
layer present between the coflow jet and the ambient air. Hence, the time scale dilemma was dealt with 
by forcing the hybrid RAS/LES model to remain in RAS mode for the outer portion of the coflow jet and 
ambient region. In particular, an initial solution was obtained for the entire flowfield using a steady-state 
RAS. The calculation was then restarted with the hybrid RAS/LES approach with the blending function 
forced to unity for r/D > 2. This value was found to be large enough to contain all convecting coherent 
structures formed in the mixing layer between the two jet flows for the LES domain of interest. 

The form of SGS closure used in this effort was based on the one-equation Yoshizawa 20 model: 
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dt 
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However, two noteworthy modifications to the model were made for this effort. The first modification 
involves the definition of the filter width, A. Traditional LES practitioners have tended to work with truly 
isotropic grids, justifying the use of a simple isotropic definition for the filter width (e.g. the cubed root of 
the cell volume). Grids utilized for configurations of engineering interest, however, will have a some level of 
anisotropy. Hence, hybrid RAS/LES practitioners tend to prefer the following anisotropic definition for the 
filter width: 

A = max (A u Aj,A k ) (7) 

In this expression, A *, Aj, and A& denote the width along each computational direction of a three- 
dimensional hexahedral grid cell. The second modification affects the level of SGS viscosity in LES regions. 
At equilibrium (defined when production balances dissipation in Eq. 6a), the Yoshizawa model reduces to an 
algebraic Smagorinsky closure model. The effective Smagorinsky constant (C s ) implied by the model under 
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Figure 9: Comparison of argon mass fraction (RAS predictions) with measured values 
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Figure 10: Comparison of Case 2 Pitot pressure (RAS predictions) with measured values 
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equilibrium conditions is given by the following expression: 

{C s f = V2C, (A) 2 (8) 

The default values given in Ref. 20 for C M and Cd are 0.05 and 1.0, respectively. The value of C s implied 
by these values is 0.126. There is no universal value for this Smagorinsky constant, but the LES community 
has typically endorsed values on the order of 0.2 for homogeneous turbulence and 0.065 for shear flows. 24 
Initial calculations used a reduced value of 0.02075 for C M to recover the “accepted” Smagorinsky constant 
for shear flows. This value, however, proved to be too dissipative within the present second-order upwind 
numerical framework. The numerical dissipation is on the same order as the SGS viscosity (A 2 ), hence one 
might expect lower values to be required. A value of C M = 0.00823 was found to readily unlock significant 
levels of resolved turbulence energy in the present simulations. This value enforces a Smagorinsky constant 
of 0.0325 (half the recommended value for shear flows) under equilibrium conditions, and was selected as 
the baseline setting. 



Figure 11: Instantaneous Schlieren and normalized helium mass fraction contours (hybrid RAS/LES) 

An instantaneous snapshot of the flowfield for Case 1 is shown in Fig. 11. The top image displays a 
numerical Schlieren of the flow structure, and the bottom image shows the normalized helium mass fraction 
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contours. A significant level of resolved turbulent content was captured in the simulation. The recirculation 
zone behind the base of the center-body provided a source of large scale unsteadiness that rapidly triggered 
Kelvin Helmholtz instabilities in the shear layer. These instabilities transitioned the flow from a fully modeled 
Reynolds- averaged state exiting the nozzle, to a primarily resolved turbulent state a few jet diameters 
downstream of the nozzle exit plane. Unsteady shock and expansion waves are seen reflecting through the 
jet structure. The shear layer between the outer jet and ambient air maintained the desired steady-state 
Reynolds- averaged behavior due to the forced usage of the RAS equations outside of r/D = 2.0. The high 
level of shear associated with this operating condition provided a strong mechanism for turbulence production 
as evident by the rich array of vortical structures seen in the helium contours. These structures provide the 
turbulent “stirring” that leads to enhanced mixing by stretching the interface between the air and injectant. 

Prior to gathering any statistics, time integrations were performed for two complete flow-through times to 
flush out the initial transients during the transition from the converged RAS solution to the hybrid RAS/LES 
result. This relatively small time interval was sufficient due to a lack of any large regions of separated flow, and 
the fact that the ambient region was essentially frozen in its Reynolds- averaged state. At this point, running 
averages were computed and stored after each iteration. In order to reduce the integration time required 
to gather meaningful statistics, a combination of temporal and spatial (in the circumferential direction) 
averaging was performed. The spatial-average was made possible by a three-step process. The first step 
was to interpolate the ensemble- averaged data onto a polar grid topology (in the cross-flow plane) at the 
experimental streamwise stations. The second step involved a coordinate transformation of the flow data to 
a cylindrical coordinate system (x,y,z — > x,r, 6). The final step was to average the transformed variables 
along the azimuthal curves of constant radius (which corresponds to the spatially homogeneous direction). 
All three steps were performed using the TECPLOT post-processing package 25 and automated using the 
TECPLOT macro scripting language. 

The averaged hybrid RAS/LES helium mass fraction profiles are compared with the measurements in 
Fig. 12. In general, the predictions show a shear layer growth rate that is more rapid than that implied 
by the measurements. As mentioned previously, the SGS model coefficients were adjusted to promote 
the onset of flow instabilities. It appears that this adjustment leads to an inappropriate level of modeled 
viscosity (i.e. a value that is too low) once the Kelvin Helmholtz instabilities have transitioned to a fully 
turbulent state. This situation might be improved if resolved turbulent structures were allowed to develop 
at the nozzle exit plane. This would remove the delay in transition that is associated with the conversion 
of modeled Reynolds-averaged turbulence to resolved LES content, and would allow a larger Smagorinsky 
coefficient to be utilized. One reasonably economic approach to this task is to implement recycling/rescaling 
techniques. 26, 27, 28 This will be the focus of future efforts. A comparison of the averaged hybrid RAS/LES 
results with those obtained using pure RAS at the default (Sc t = 0.5) condition show many similarities. This 
particular RAS simulation also over-predicted the shear layer growth rate, but one noteworthy difference 
between the results is the manner in which the helium mass fraction asymptotically approaches zero at the 
edge of the mixing layer. RAS approaches are notorious for producing a sharp non-physical interface at the 
edge of the shear layer. The hybrid RAS/LES results show a much smoother asymptotic decay that is more 
representative of the measurements. Statistical convergence has been assessed by comparing the statistics 
extracted over successively larger iteration intervals. The statistics extracted over twenty and thirty thousand 
iterations show practically identical results. Another useful statistical convergence check is to examine the 
averaged azimuthal velocity component. This value should approach zero as the statistics are converged. 
The maximum absolute value of azimuthal velocity (averaged over 30,000 cycles) was 2.5 m/s, or roughly 
0.2% of the streamwise velocity. 

The averaged Pitot pressure profiles are compared with the measurements in Fig. 13. These results 
reinforce the assertion that the mixing layer growth rate has been over-predicted, particularly at stations 
downstream of an x/T of 4. Although the deficit in Pitot pressure profile is broader across the shear layer, 
the qualitative features of the profile shape mimics that of the measurements. This was not the case in the 
RAS results that over-predicted the level of mixing. The finer details of the flow structure, e.g. the slight 
pressure deficit due the the Mach disk at the axis, are also evident in the hybrid RAS/LES data. As was the 
case with the helium mass fraction statistics, the averaged Pitot pressure profiles for both averaging intervals 
show practically identical results, suggesting that the first-order moments are well converged. 

An instantaneous snapshot of the flowfield for Case 2 is shown in Fig. 14. The top image displays a 
numerical Schlieren of the flow structure, while the lower image displays the argon mass fraction distribution. 
The Kelvin Helmholtz instabilities are clearly evident in both images and persist for many jet diameters before 
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Figure 12: Comparison of normalized helium mass fraction (averaged hybrid RAS/LES predictions) with 

measured values 
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Figure 13: Comparison of Case 1 Pitot pressure (averaged hybrid RAS/LES predictions) with measured values 
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breaking down into fully three-dimensional turbulent structures. This transition process occurs near an x/D 
of 10 or 12, where the last visible wave structure appears in the Schlieren image. The low level of shear for 
this case (as compared with Case 1) is the root cause for the delay in the break-down of these structures. 
Large (i.e. on the order of the jet diameter) turbulent coherent structures do not form until an x/D of 15 to 
20. Hence, the turbulent “stirring” processes that lead to enhanced mixing are not prevalent until the latter 
stages of the LES domain. 

The averaged hybrid RAS/LES argon mass fraction profiles are compared with measurements in Fig. 15. 
Overall, the level of turbulent mixing between the argon and air jets was under-predicted. These results 
are comparable with those obtained using the RAS approach (see Fig. 9), though the hybrid RAS/LES 
shows a slightly reduced mixing level at the last experimental station. The smooth asymptotic transition 
of the argon profile towards zero at the edge of the mixing layer is again evident in the hybrid RAS/LES 
results. The reduced spreading rate of the shear layer is also seen in the Pitot pressure profile comparisons 
displayed in Fig. 16. The peak shear layer deficit is consistently over-predicted at each station. Future 
work will examine the influence of providing realistic resolved turbulent content exiting the coaxial nozzles. 
This may be particularly important for the present application because the boundary layer thickness of the 




Figure 14: Instantaneous Schlieren and argon mass fraction contours (hybrid RAS/LES) 

approach coflow is an order of magnitude thicker than the blunt base region which is currently triggering 
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the unsteadiness. The statistics gathered after 10,000 cycles show practically identical results as the data 
averaged over 20,000 cycles, suggesting that the statistics are converged for these first-order correlations. 
The maximum azimuthal velocity (averaged over 20,000 cycles) was 0.6 m/s. 

Second-Order Correlations 

The primary second-order correlations that are modeled in the RAS equation set are the Reynolds stress 
tensor and the Reynolds heat and mass flux vectors. The Reynolds stress tensor is typically modeled by the 
Boussinesq approximation: 


pu/u- 



pk + fit 


du k \ 

dx k ) 




and the Reynolds flux vectors commonly utilize gradient diffusion models, i.e. 
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The accuracy of the modeling employed for these terms is critical for the success of RAS approaches. Mea- 
sured data for each of the above correlations is scarce, so it is often difficult to directly assess the accuracy 
of the RAS closures. One of the primary goals of this effort was to use the hybrid RAS/LES simulations 
to assess the performance of the RAS modeling. Ideally, one would hope for a close match between the 
hybrid RAS/LES results and available measurements, but the level of agreement obtained in this effort was 
not entirely satisfying. Nevertheless, some insight into the RAS modeling can be gleaned by comparing the 
modeled terms with values extracted from the resolved LES field. 

The rms of the streamwise velocity fluctuation is considered first. This quantity was the only second- 
order correlation that was measured, and the measurement was only taken for Case 1. These measurements 
were made using the RELIEF 12 oxygen flow-tagging technique. Figure 17 compares the measured values 
with the values computed from the RAS model (Eq. 9, Sc t = 0.5) as well as the ensemble-averages extracted 
from the hybrid RAS/LES. The modeled RAS results agree remarkably well with the measurements in 
the early stages of the shear layer development, but the width of the profile is under-predicted at stations 
further downstream. The hybrid RAS/LES results over-predict the rms values in the shear layer at the first 
station, but the profiles at the stations further downstream compare favorably with the measurements. One 
particular noteworthy feature is the build-up of turbulent fluctuations in the core flow near the axis that 
is evident in the hybrid RAS/LES results. This feature is absent from the RAS result because the mean 
flow gradients are relatively small outside of the shear layer. Mean velocity gradients are the sole source of 
turbulence production for most RAS models. The enforcement of a RAS behavior for r/D > 2 in the hybrid 
RAS/LES simulations prevented the formation of unsteady turbulence structures in the outer jet/ambient 
air shear layer; limiting the build-up of turbulence in the core of the outer jet. Hence the core flow values in 
the outer jet, while larger than that predicted by pure RAS, is under-predicted. 

The rms radial and azimuthal velocity fluctuations are compared in Figs. 18 and 19. Linear eddy viscosity 
models based on the Boussinesq approximation tend to produce nearly isotropic velocity variances (normal 
stresses) when applied to strain dominated flows. The degree of anisotropy is governed by the magnitude of 
the velocity gradients (other than those contributing to the strain rate) relative to 2/3 pk term (see Eq. 9). 
Hence, it is not surprising that the RAS model returned a nearly isotropic set of normal stress tensors for 
this turbulent shear flow. The hybrid RAS/LES data, on the other hand, showed a non- negligible level of 
anisotropy in the velocity variances. In particular, the streamwise velocity variance was substantially larger 
than the cross-flow variances. This result is typical for shear dominated flows. 

The Reynolds shear stress ( u"v ") and Reynolds mass flux Y^ e v" terms are compared in Figs. 20 and 21. 
The accuracy of the RAS equations, to a large extent, is driven by how well these terms are modeled. Since 
the Sc t = 1.0 RAS results compared favorably to the measured mean flow properties, correlations from this 
simulation were also added to the figures. The shear stress extracted from the hybrid RAS/LES data is larger 
than those produced by the RAS models at the first axial two stations. The situation is reversed at the last 


20 of 27 


American Institute of Aeronautics and Astronautics 



2.0 


a Measured 

1.5 


- RAS/LES (20000 cycles) 
RAS/LES (10000 cycles) 

9 i.o 



n c; 

L* 

x/D = 0.9969 

u.o 


1 ... 1 ... 1 ... 1 . i < 1 


,qI 1 1 1 * 

0.0 0.2 0.4 0.6 0.8 1.0 









Figure 15: Comparison of argon mass fraction (averaged hybrid RAS/LES predictions) with measured values 
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Figure 16: Comparison of Case 2 Pitot pressure (averaged hybrid RAS/LES predictions) with measured values 
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two stations. This result is consistent with the fact that the hybrid RAS/LES results predicted a more rapid 
mixing rate than the RAS approaches. A similar behavior is evident in the Reynolds mass flux vector. The 
hybrid RAS/LES data showed the most rapid mixing, the Sc t = 1.0 RAS yielded the slowest mixing rate, 
and the Sc t = 0.5 RAS result was somewhere inbetween. Hence at the early axial stations, the turbulent 
mixing term should be largest for the hybrid RAS/LES and smallest for the Sc t = 1.0 RAS. Eventually, the 
condition with the most efficient mixing will show a decayed mixing rate since a large fraction of the mixing 
process has completed. Thus at some station downstream, the simulation with the slower mixing rate should 
exceed that of the more rapid mixing rate. This is precisely the situation that occurs between the second 
and third streamwise stations shown in Fig. 21. 






Figure IT: Comparison of Case 1 streamwise velocity fluctuation rms with measured values 
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Figure 18: Comparison of Case 1 radial velocity fluctuation rms 






Figure 19: Comparison of Case 1 azimuthal velocity fluctuation rms 
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Figure 20: Comparison of Case 1 u v" (Reynolds shear stress correlation) 






Figure 21: Comparison of Y^ e u” (Reynolds mass flux correlation) 
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Summary 


Reynolds- averaged and hybrid Reynolds- averaged/large-eddy simulations have been performed to model 
a supersonic coaxial jet flow experiment. The experiment consists of an outer jet of air, and an inner jet 
that was either a He — O 2 mixture or pure argon. Both jets exhausted into an ambient environment. The 
Mach 1.8 nozzle flows exiting the test apparatus were nearly pressure- matched for both injectant conditions. 
However, the level of shear and the compressibility of the mixing layer varied depending on the injectants 
involved. The helium condition resulted in a highly compressible mixing layer with a convective Mach 
number of 0.7. The argon condition, on the other hand, was nearly velocity- matched and produced a mixing 
layer with a convective Mach number of 0.16. A comprehensive set of measurements were taken which 
include: Pitot pressure, mean and rms velocities, and gas sampling. The model geometry, flow conditions, 
and measurement uncertainties were all well documented, resulting in a package that is well suited for model 
validation efforts. 

The goal of the computational effort was to assess the state-of-the-art for both RAS and hybrid RAS /LES 
approaches as applied to compressible turbulent flows. The Reynolds- averaged simulation for the helium 
cases displayed a strong sensitivity to choice of turbulent Schmidt number. A value of 1.0 was found to be 
an optimal choice for this flow condition. A lower turbulent Schmidt number, on the order of 0.5, provided 
the best match with measurements for the argon case, however. The uncertainty involved with appropriate 
choices for RAS modeling parameters highlights the difficulty with using these approaches for predictive 
simulations. In principal, LES or hybrid RAS/LES methods have the potential to reduce this uncertainty 
by resolving a substantial fraction of the turbulent flowfield. The hybrid simulations performed in this 
effort, however, were no more predictive than the baseline Reynolds- averaged predictions. The explanation 
provided for the discrepancy between the hybrid RAS/LES results and the measurements centered around 
issues related to how the Reynolds- averaged state transitions to a state with resolved turbulent content. The 
addition of resolved turbulent content to the inflow conditions was suggested to address this concern. Finally, 
comparisons made between resolved second-order turbulence statistics and their modeled Reynolds-averaged 
counterparts were discussed. Refined hybrid simulations are required, however, to improve the accuracy of 
the hybrid RAS/LES results before any solid conclusions can be drawn concerning the RAS modeling. 
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